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AN UNCONDITIONNALLY STABLE PRESSURE CORRECTION SCHEME FOR 
COMPRESSIBLE BAROTROPIC NAVIER-STOKES EQUATIONS 

T. Gallouet^, L. Gastaldo^, R. Herein^ and J.C. Latche^ 

Abstract. We present in this paper a pressure correction scheme for barotropic compressible Navier- 
Stokes equations, which enjoys an unconditional stability property, in the sense that the energy and 
maximum-principle-based a priori estimates of the continuous problem also hold for the discrete solu- 
tion. The stability proof is based on two independent results for general finite volume discretizations, 
both interesting for their own sake: the L'^-stability of the discrete advection operator provided it is 
consistent, in some sense, with the mass balance and the estimate of the pressure work by means of 
the time derivative of the elastic potential. The proposed scheme is built in order to match these 
theoretical results, and combines a fractional-step time discretization of pressure-correction type to 
a space discretization associating low order non-conforming mixed finite elements and finite volumes. 
Numerical tests with an exact smooth solution show the convergence of the scheme. 
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1. Introduction 

The problem addressed in this paper is the system of the so-called barotropic compressible Navier-Stokes 
equations, which reads: 

|f+V.(p.) = 

^(pw) + V • Vp- V • t(m) = / ^^'^ 

p q{p) 

where t stands for the time, p, u and p are the density, velocity and pressure in the flow, / is a forcing term and 
t{u) stands for the shear stress tensor. The function g{-) is the equation of state used for the modelling of the 
particular flow at hand, which may be the actual equation of state of the fluid or may result from assumptions 
concerning the flow; typically, laws as p = P^^'^, where 7 is a coefficient specific to the fluid considered, are 
obtained by making the assumption that the flow is isentropic. This system of equations is posed over 57 x (0, T), 
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where O is a domain of W^, d < 3 supposed to be polygonal (d = 2) or polyhedral {d = 3), and the final time T 
is finite. It must be supplemented by boundary conditions and by an initial condition for p and u. 

The development of pressure correction techniques for compressible Navier-Stokes equations may be traced 
back to the seminal work of Harlow and Amsden [20, 21] in the late sixties, who developped an iterative 
algorithm (the so-called ICE method) including an elliptic corrector step for the pressure. Later on, pressure 
correction equations appeared in numerical schemes proposed by several researchers, essentially in the finite- 
volume framework, using either a collocated [10, 24, 26, 30, 33, 34] or a staggered arrangement [2, 4, 7, 22, 23, 
25,37,38,40-42] of unknowns; in the first case, some corrective actions are to be foreseen to avoid the usual 
odd-even decoupling of the pressure in the low Mach number regime. Some of these algorithms are essentially 
implicit, since the final stage of a time step involves the unknown at the end-of-step time level; the end-of- 
step solution is then obtained by SIMPLE-like iterative processes [10,24-26,30,34,39]. The other schemes 
[2,7,22,23,33,37,38,40,42,43] are predictor-corrector methods, where basically two steps are performed in 
sequence: first a semi-explicit decoupled prediction of the momentum or velocity (and possibly energy, for non- 
barotropic flows) and, second, a correction step where the end-of step pressure is evaluated and the momentum 
and velocity are corrected, as in projection methods for incompressible flows (see [5,36] for the original papers, 
[29] for a comprehensive introduction and [19] for a review of most variants). The Characteristic-Based Split 
(CBS) scheme (see [31] for a recent review or [44] for the seminal paper), developped in the finite-element 
context, belongs to this latter class of methods. 

Our aim in this paper is to propose and study a simple non-iterative pressure correction scheme for the 
solution of ll]). In addition, this method is designed so as to be stable in the low Mach number limit, since 
our final goal is to apply it to simulate by a drift-flux approach a class of bubbly flows encountered in nuclear 
safety studies, where pure liquid (incompressible) and pure gaseous (compressible) zones may coexist. To this 
purpose, we use a low order mixed finite element approximation, which meets the two following requirements: 
to allow a natural discretization of the viscous terms and to provide a spatial discretization that is intrinsically 
stable {i.e. without the adjunction of stabilization terms to circumvent the so-called inf-sup or BB condition) 
in the incompressible limit. 

In this work, a special attention is payed to stability issues. To be more specific, let us recall the a priori 
estimates associated to problem ([Ij with a zero forcing term, i.e. estimates which should be satisfied by any 
possible regular solution [15,28,32]: 

{i) p{x,t) > 0, Vx- e 17, e (0,r) 



(2) 



(m) / p(a;,t)da;== / p(a;, 0) dx, Vte(0,r) 
Jn Jn 

\ d f d f 

{Hi) 2 di J Pl^;,*)""!-^,*)^ da; + — J p{x,t) P{p{x,t)) dx 

+ / t{u{x, t)) : \'u{x, t) d.T = 0, Mte (0, T) 

In the latter relation, P, the "elastic potential", is a function derived from the equation of state, which satisfies: 

P'{z) = ^ (3) 

where p{-) is the inverse function of p(-), i.e. the function giving the pressure as a function of the density. The 
usual choice is, provided that this expression makes sense: 

p(z) =r^as (4) 
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For these estimates to hold, the condition ©-(i) must be satisfied by the initial condition; note that a non-zero 
forcing term / in the momentum balance would make an additional term appear at the right hand side of 
relation ^-{iii). This latter estimate is obtained from the second relation of ([T|) {i.e. the momentum balance) 
by taking the scalar product by u and integrating over fl. This computation then involves two main arguments 
which read: 



d 

— {pu) + \/-{pu®u) 



• u da: = - — / pu^ dx 



(z) Stability of the advection operator: 

(a) Stability due to the pressure work: [ —p V ■ udx ^ — [ p(x, t) P{p{x, t)) dx 

Jn at Jo 



(5) 



Note that the derivation of both relations make use of the fact that the mass balance equation holds in a crucial 
way. 

This paper is organized as follows. 

As a first step, we derive a bound similar to ©-(m) for a given class of spatial discretizations; the latter are 
introduced in section 12.11 and the desired stability estimate (theorem 12. ip is stated and proven in section 12.21 
We then show that this result allows to prove the existence of a solution for a fairly general class of discrete 
compressible flow problems. Section [5] gathers this whole study, and constitutes the first part of this paper. 

In a second part (section [3]) , we turn to the derivation of a pressure correction scheme the solution of which 
satisfies the whole set of a priori estimates To this purpose, besides theorem 12.11 we need as a second key 
ingredient a discrete version of the bound ©-(i) relative to the stability of the advection operator, which is 
stated and proven in section [221 ftheorem l3.ip . We then derive a fully discrete algorithm which is designed to 
meet the assumptions of these theoretical results, and establish its stability. Some numerical experiments show 
that this scheme seems in addition to present, when the solution is smooth, the convergence properties which 
can be expected, namely first order in time convergence for all the variables and second order in space in 
and discrete norm for the velocity and the pressure, respectively. 



2. Analysis of a class of discrete problems 



The class of problems addressed in this section can be seen as the class of discrete systems obtained by space 
discretization by low-order non-conforming finite elements of continuous problems of the following form: 

Au + Vp ^ f in 

^H^+V.(,(p).)=0 inn (6) 
u = on dfl 

where the forcing term / and the density field p* are known quantities, and A stands for an abstract elliptic 
operator. The unknown of the problem are the velocity u and the pressure p; the function g(-) stands for the 
equation of state. The domain is a polygonal (d = 2) or polyhedral {d = 3) open, bounded and connected 
subset of M''. Of course, at the continuous level, this statement of the problem should be completed by a 
precise definition of the functional spaces in which the velocity and the pressure are searched for, together with 
regularity assumptions on the data. This is out of context in this section, as system ^ is given only to fix 
ideas and we restrict ourselves here to prove some mathematical properties of the discrete problem, namely to 
establish some a priori estimates for its solution and to prove that this nonlinear problem admits some solution 
for fairly general equations of state. 

This section is organized as follows. We begin by describing the considered discretization and precisely stating 
the discrete problem at hand. Then we prove, for the particular discretization at hand, a fundamental result 
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which is a discrete analogue of the elastic potential identity. The next section is devoted to the proof of the 
existence of a solution, and we finally conclude by giving some practical examples of application of the abstract 
theory dcvcloppcd here. 

2.1. The discrete problem 

Let be a decomposition of the domain il cither into convex quadrilaterals {d — 2) or hexahedrons {d = 3) 
or in simplices. By £ and £{K) we denote the set of all {d — l)-edges a of the mesh and of the element K ^ A4 
respectively. The set of edges included in the boundary of ft is denoted by fext and the set of internal ones 
(i.e. £ \ £cxt) is denoted by £int- The decomposition M is supposed to be regular in the usual sense of the 
finite element literature (e.g. [6]), and, in particular, M satisfies the following properties: (} = [Jkem^' 
K, L G A4, then _R'nZ = 0ori?nLisa common face of K and L, which is denoted by K\L. For each internal 
edge of the mesh a — K\L, ukl stands for the normal vector of a, oriented from K to L. By \K\ and \a\ we 
denote the measure, respectively, of K and of the edge a. 

The spatial discretization relics cither on the so-called "rotated bilinear element" /Pq introduced by Rannacher 
and Turek [35] for quadrilateral of hexahedric meshes, or on the Crouzeix-Raviart element (see [8] for the seminal 
paper and, for instance, [12, p. 83-85] for a synthetic presentation) for simplicial meshes. The reference element 
K for the rotated bilinear element is the unit c?-cube (with edges parallel to the coordinate axes); the discrete 
functional space on K is Qi{KY, where Qi{K) is defined as follows: 

(3i(ii:) = spanjl, {xi)i=i„„^d, {xf - x^^^)t=i d-i} 

The reference element for the Crouzeix-Raviart is the unit d-simplex and the discrete functional space is the 
space Pi of affine polynomials. For both velocity elements used here, the degrees of freedom are determined by 
the following set of nodal functionals: 

{F„.„ ae£{K),i = l,...,d}, F„,,{v) = \n\-^ f (7) 

The mapping from the reference clement to the actual one is, for the Rannacher- Turck clement, the standard 
Q\ mapping and, for the Crouzeix-Raviart clement, the standard affine mapping. Finally, in both cases, the 
continuity of the average value of discrete velocities (i.e., for a discrete velocity field u, Fcr,i(w)) across each face 
of the mesh is required, thus the discrete space Wh is defined as follows: 

= { w/i € L^(ri) : Vh\K e Q\{KY , \/K G M] Fc^i{vh) continuous across each edge a e £int, 1 < « < rf; 

F„^,{vh) = 0, Vd G fext, 1 < « < d } 

For both Rannacher- Turck and Crouzeix-Raviart discretizations, the pressure is approximated by the space Lh 
of piecewise constant functions: 

Lh = {qh e L^{^) : qh\K ^ constant, V/i G 7W} 

Since only the continuity of the integral over each edge of the mesh is imposed, the velocities are discontinuous 
through each edge; the discretization is thus nonconforming in H^{^lY. These pairs of approximation spaces 
for the velocity and the pressure are inf-sup stable, in the usual sense for "piecewise H^" discrete velocities, i.e. 
there exists Ci > independent of the mesh such that: 

j p V • da; 

Vp G Lh, sup "'^11 11 > Ci \\p - p||L2(n) 

v(^Wh \\V\\l,b 
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where p is the mean value of p over 11, the symbol / stands for \J / and || • ||i stands for the broken 
Sobolev semi-norm: 

I |Vwpdx:= I \\/v\^dx 

l^r- hA J K Jn.h 



\l,b — 



From the definition ([7]) , each veloeity degree of freedom can be univoquely associated to an element edge. We 
take benefit of this correspondence by using hereafter, somewhat improperly, the expression "velocity on the 
edge ct" to name the velocity vector defined by the degrees of freedom of the velocity components associated to 
a. In addition, the velocity degrees of freedom are indexed by the number of the component and the associated 
edge, thus the set of velocity degrees of freedom reads: 

{Vcr,i, a e £int, i <i < d} 

We define = J2i=i where a is the i*'' vector of the canonical basis of M''. We denote by ip^^ the vector 

shape function associated to Va,i, which, by the definition of the considered finite elements, reads: 

(i) 

where is a scalar function. Similarly, each degree of freedom for the pressure is associated to a mesh K, and 
the set of pressure degrees of freedom is denoted by {px, K G A4}. 

For any K G M, let p*j^ be a quantity approximating a known density p* on K. The family of real numbers 
{p*k)k^m is supposed to be positive. The discrete problem considered in this section reads: 

a(u, (^W) - / p V • (/pW = / / • ^« Ax Va G £int {<J = K\L), 1 < ^ < d 

Jnji Jn 

\K\ ^^P"^^-P*^' + J2 <K b{pk) v-^ g{pL) ^0 VKeM ^ 

(j=K\L 

where vj^^ and v^^ stands respectively for vj^^ = max(vCT,ii:, 0) and v^^ — max(— Vo.,;^, 0) with Va.K = 
\(j\ucr ■ riKL " K ^ '^a K- fi^st equation is the standard finite element discretization of the first relation 
of dS]), provided the following identity holds: 

Vn G Wh, Vui G Wh, a{v, w) — Av ■ wdx 

Jn 

As the pressure is piecewisc constant, the finite element discretization of the second relation of i.e. the mass 
balance, is similar to a finite volume formulation, in which we introduce the standard first-order upwinding. 
The bilinear form a(-, •) is supposed to be elliptic on Wh, i.e. to be such that the following property holds: 



3ca > such that, Vw G Wh, a{u, u) > Ca ||u|| 



2 



where || • ||* is a norm over Wh. We denote by || • ||* its dual norm with respect to the L'^{Q)'^ inner product, 
defined by: 



V . wdx 

n 



G Wh, \\v\\* = sup ,, ,, 
wew,, \\w\\ 
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2.2. On the pressure control induced by the pressure forces work 

The aim of this subsection is to prove that the discretization at hand satisfies a stabihty bound which can 
be seen as the discrete analogue of equation ©-(m), which we recall here: 

I pV ■udx= ^ I pP{p)dx, where P(-) is such that P'{z) = 

The formal computation which allows to derive this estimate in the continuous setting is the following. The 
starting point is the mass balance, which is multiplied by [pP{p)y: 



This relation yields: 



And thus: 



d[pp{p)] 



= 



dt 

d[pPip)] 



Developping the derivative, we get: 

d[pP{p)] 
dt 



dt 



+ [pPip)]' [u • Vp + pV • u] = 
u ■ V[pP(p)] + [pPip)]'pV • u = 



(9) 



V[pP(p)] + pP(p)V • u + p2[p(p)]'v • u = 



(10) 



V-(pP(p)w) 



pV • u 



and the result follows by integration in space. We are going here to reproduce this computation at the discrete 
level. 

Theorem 2.1 (Stability due to the pressure work). Let the elastic potential P defined by ^ be such that the 
function z ^ zP[£) is once continuously differentiable and strictly convex. Let {pk)ki^m '^^'^ {Pk)ki^m satisfy 
the second relation of ([5]). Then the following estimate holds: 



f pW-udx=Yl ~P'< E ''^''<^Tt E \K\[pkP{pk)^P*kP{Pk)] 



(11) 



a=K\L 



Proof. Let us write the divergence term in the discrete mass balance over K [i.e. the second relation of ([8|) 
under the following form: 

^ PaVa,K 
a=K\L 

where is either px = q{pk) if ^a.K > or = q{Pl) if ^(j,k < 0. Multiplying this term by the derivative 
with respect to p of pP{p) computed at pK, denoted by [pP(p)]p^, we obtain: 



Pdiv,K = [PP(P)];^ PayaJ<^[pPip)]'p, 
a=K\L 



{Pa - PK)VaJ< + PK ^ V^j^ 
a=K\L a=K\L 



This relation is a discrete equivalent to equation ©: up to the multiplication by l/lifl, the first summation 
in the right hand side is the analogue of u ■ Vp and the second one to pV ■ u. Developping the derivative, we 
obtain the equivalent of relation PU)) : 

Tdiv^K =[pP{p)]'p^ Y {Pa - PK)VaJ< + PK P{Pk) Y1 ^'^:K + Pk P' {Pk) Y1 ^"^^ 
cr=K\L cr=K\L a=K\L 
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By definition ^ of P, the last term reads pk J2a=K\L ^cr,K- The process will be completed if we put the first 
two terms in divergence form. To this end, let us sum up the T(iiv,if and reorder the summation: 

TdivJ< ^ ^ Pk ^ Vcr^K + ^ ?div,o- (13) 
KGM KGM a=K\L ae£i^t 

where, if cr = K\L: 

Tdiv,<T = v^j^ [pkP(.Pk) + [pP(.p)]'p,^{Pa - pk) - PlP{pl) - [p P{p)]'p^{Pa - Pl)] 

In this relation, there are two possible choices for the orientation of cr, i.e. K\L or L\K^ where K and L are two 
cells such that cr = n L; we choose this orientation to have Va,K > 0. Let pa be defined by: 

if /OK ^ PL : PKP{pK) + [pP{p)]'p^{Pa- Pk)= PlP{pl) + [pP{p)]'pSp-- Pl) 

(14) 

otherwise : pa = pk = Pl 

As the function z ^-^ z P{z) is supposed to be once continuously differentiable and strictly convex, the technical 
lemma [^31 proven hereafter applies and pa is uniquely defined and satisfies p^ S [min(p/^, p^), max(/9/< , p^)]. 
By definition, the choice p^ = pa is such that the term Tdiv.tr vanishes, which means that, by this definition, we 
would indeed have obtained that the first two terms of equation (|12p are a conservative approximation of the 
quantity V • pP{p)u appearing in equation pO)l . with the following expression for the flux: 

Fa,K ^ [pP{p)]<jVcr,K, with: 

[pP{p)]a = PkP{pk) + [pP{p)]'p^-{Pa - pk) = PlP{pl) + [p P{p)]'pAP- - Pl) 
Whatever the choice for pa may be, we have: 

Tdiv,. - y.,K iPa - P.) {[pPip)]'p, - [pPip)]'pJ 

With the orientation taken for cr, an upwind choice yields: 

Tdiv,. = V,,K ipK - pa) i[pP{p)rp, - [pP{p)]'pJ 

and, using the fact that z i— > [pP{p)]'^ is an increasing fuction since z i-^ zP{z) is convex and that pa € 
[min(/9xj Pl), inax(p;<-, pi,)], it is easily seen that Tdiv.o- is non-negative. 

Multiplying by [pP{p)]'pj^ the mass balance over each cell K and summing for K ^ M. thus yields, by 
equation P^ : 

- E E ^'^.K = R+ E i^^i [pnp)U (15) 

K^M cr=K\L KeM 

where R is non-negative, and the result follows invoking once again the convexity of z i-^ z P{z). □ 

Remark 2.2 (On a non-dissipative scheme). The preceding proof shows that, for a scheme to conserve the 
energy {i.e. to obtain a discrete equivalent of ^-{iii)), besides other arguments, the choice of pa given by (fT4|) 
for the density at the face of the control volume in the discretization of the flux in the mass balance seems to be 
mandatory; any other choice leads to an artificial dissipation in the work of the pressure forces. Note however 
that, this discretization being essentially of centered type, the positivity of the density is not warranted in this 
case. 

In the course of the preceding proof, wc used the following technical lemma. 
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Lemma 2.3. Let g{-) be a strictly convex and once continuously derivable real function. Let a be an internal 
edge of the mesh separating the cells K and L. Then the following relations: 



ifPK^PL ■■ 9{Pk) + g' {pK){Pcr ~ Pk) = g{PL) + g' {PL){Pa - Pl) 

otherwise : ~ pK ~ PL 



(16) 



uniquely defines the real number p„. In addition, we have pn € [min(p/(-, p^), niax(pif , p^)]. 

Proof. If pk — PL, there is nothing to prove. In the contrary case, without loss of gcnerahty, let us choose K 
and L in such a way that pK < Pl- By reordering equation (|16p . we get: 

9{pk) + g'{pK) {.PL - pk) - g{pL) = {pa - pl) [g'{pL) ~ g'{pK)] 

As, because (?(•) is strictly convex, g' {pl) — g'ipx) does not vanish, this equation proves that pa- is uniquely 
defined. In addition, for the same reason, the left hand side of this relation is negative and g' {pl) — g' {pk) is 
positive, thus we have Pa < Pl- Reordering in another way equation (|16p yields: 



g{pL) + 9'{pl) (pk - pl) - 9{pk) = {pa - pk) [g'ipK) - g'ipL)] 

which, considering the signs of the left hand side and of g'{pK) ^ g'iPL), in turn implies p^ > px- 
2.3. Existence of a solution 



□ 



The aim of this section is to prove the existence of a solution to the discrete problem under study. It follows 
from a topological degree argument, linking by a homotopy the problem at hand to a linear system. 

This section begins by a lemma which is used further to prove the strict positivity of the pressure. 
Lemma 2.4. Let us consider the following problem: 



St 



^a,K'P2{pK)-y^^K'P2{pL) =0 



(17) 



(T=K\L 



where ipi is an increasing function and ipi is a non- decreasing and non-negative function. Suppose that there 
exists p such that: 



(pi{p) + St(p2{p) max 
KeM 



Va.K 



<j=K\L 



= Tmu\Lpi{p\)] 
KeM 



(18) 



Then, \/K G M, pK satisfies pk > P- 

Proof. Let us assume that there exists a cell K such that pf^ = miiiKeMPK < P- Multiplying by (5t/|A'| the 
relation ([TT]) written for K = K, we get: 



St 



bi?) + 7^ (v^,K^2(Pi?)-V^_^^2(pL)) =¥'l(pjf) 



(19) 



a=K\L 



Then substracting from (|18p we have: 



St 



"PiiPK) ~ V'i(p) + 7^ Y V2{pk) 



.K 



^2{Pl) 



=K\L 



-Stifioip) max 

K£M 



ir\ ^ 



K 



Vi(Pk) - min [</'i(Pk)] > 



KeM 
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The previous relation can be written as Ti +T2 +T3 > with: 

Ti = fiiPR) - flip) 



T2 = 6tip2{pK) 

6t 



y ^ 



^K\L 



St ifio (p) max 



0' E ^--^.^ 

cT=K\L ' ' 



T: 



\K\ 



a=K\L 



As ipi{-) is an increasing function and, by assumption, pj^ < p, we have Ti < 0. Similarly, < if2{PK) < V2(p) 
and the discrete divergence over K [i.e. ^/\K\ X]cr=i?|L is necessarily smaller than the maximum of this 
quantity over the cells of the mesh, thus T2 < 0. Finally, as, by assumption, pf^ < p^ for any neighbouring 
cell L of K, (p2{-) is a non-decreasing function and v~^ > 0, T3 < 0. We thus obtain a contradiction with the 
initial hypothesis, which proves px > P, e A^. □ 

We now state the abstract theorem which will be used hereafter; this result follows from standard arguments 
of the topological degree theory (see [9] for an exposition of the theory and [13] for another utilisation for the 
same objective as here, namely the proof of existence of a solution to a numerical scheme). 

Theorem 2.5 (A result from the topological degree theory). Let N and M be two positive integers and V be 
defined as follows: 



V = {{x,y)e 



such that y > 0} 



where, for any real number c, the notation y > c means that each component of y is greater than c. Let 



pN 



pM 



and f and F be two continuous functions respectively from V and V x [0, 1] to 



b e 

satisfying: 

(i) F(., !) = /(.); 

(ii) Va € [0, 1], if V is such that F{v, a) ^ b then v G W , where W is defined as follows: 

W = {(x, y) e X K^^ s.t. \\x\\ < Ci and e < y < C2} 

with Ci, C2 and e three positive constants and \\ ■ \\ is a norm defined over R^; 

(iii) the topological degree of F(-,0) with respect to b and W is equal to do 7^ 0. 

Then the topological degree of F{-, 1) with respect to b and W is also equal to do 7^ 0; consequently, there exists 
at least a solution v W such that f{v) = b. 

We are now in position to prove the existence of a solution to the considered discrete Darcy system, for fairly 
general equations of states. 

Theorem 2.6 (Existence of a solution). Let us suppose that the equation of state q{-) is such that: 

(1) q{-) is increasing, q{Qi) ~ and lim g{z) = +00, 

z — '+00 

(2) there exists an elastic potential P such that (0) holds, the function z 1-^ z P{z) is once continuously dif- 
ferentiate and strictly convex and zP{z) > ~Cp, Vz G (0,-|-c»), where Cp is a non-negative constant. 

Then, there exists a solution {uc,PK)a£Eir,t,K£M to the discrete problem at hand 



Proof. This proof is obtained by applying theorem 12.51 Let N — d card(£int), M — card(AI) and the mapping 



F: V x[0,l] 



pN 



pM 



be given by: 



F{u,p, a) 



a e £int, I <i < d 



qK, K e M 
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with: 

v„^^ = a{u, .^W) _ Q, / p V • (/jW - / / • (^W 



-'n.h - , , 

\K\ _ (20) 

(j=K\L 

where, G M, p*j^ is chosen such that p*j^ = q(j)*x)', note that, by assumption, g{-) is one to one from (0, +oo) 
to (0, +oo), so tire preceding definition make sense. The problem F{u,p, 1) = is exactly the system ([5]). 

The present proof is built as follows: by applying theorem l2.H we obtain a control on u in the discrete norm, 
uniform with respect to a. Since we work on a finite dimensional space, we then obtain a control on p in L°° 
by using the conservativity of the system of equations. For the same reason, the control on u yields a bound in 
L°° of the value of the discrete divergence, which is shown to allow, by lemma [231 to bound p away from zero 
independently of a. The proof finally ends by examining the properties of the system F{u,p, 0) = 0. 

Step 1 : a E [0,1], \\ ■ ||* estimate for the velocity. 

Setting Ug. = in (j20p . multiplying the corresponding equation by u^^i and summing over a € Sint and 1 < i < d 
yields the following equation: 

a{u,u)~a / pV-itdx= / f-udx 



inji 

By a computation very similar to the proof of theorem 12. 11 we see that, from the second relation of (PO)) with 
to = 0: 



Jnji ot ,. 



where pK = Q{Pk)- By the stability of the bilinear form a(-, •) and Young's inequality, we thus get: 

y + i E 1^1 pi< p^p'<) <^\\fr^ +IiT.\^\pk np*K) (21) 

^ ' KeM KeM 



T2 

By assumption, T2 > — Cp|r2| and we thus get the following estimate on the discrete norm of the velocity: 

\\u\U < Ci (22) 
where Ci only depends on the data of the problem, i.e. 5t, f and p* and not on a. 
Step 2: a E (0, 1], L°° estimate for the pressure. 

Let us now turn to the estimate of the pressure. By conservativity of the discrete mass balance, it is easily seen 
that: 

KeM KeM 
As each term in the sum on the left hand side is non-negative, we thus have: 

which, as, by assumption, hm g{z) = +00, yields: 

z^-\-oo 



yK eM, PK < C2 (23) 
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where C2 only depends on the data of the problem. 

Step 3: a Cz (0, 1], p bounded away from zero. 
Applying lemma [513] with fii') = y^2i') ~ we get: 

yKeM, PK>Pa 

where pa is given by: 

QiPa) = 



min p(p*) 
KeM 



1 + 6t max 
KeM 



E w\ 



<y=K\L 

Note that Pa is well defined since, by assumption, q[-) is one to one from (0, +00) to (0, +00). As a < 1, we 
get: 

min p(p*) 
(- \ KeM^^"^ 



1 + St max 



0' E ^ 



=K\L 



and, by equivalence of norms in a finite dimensional space, the bound (j22p also yields a bound in the h°° norm 
and, finally, an upper bound for the right hand side of this relation. We thus get, still since q{-) is increasing 
on (0, +00) that, Va G (0, 1], Pq > ei, and, finally: 



where e only depends on the data. 
Step 4 ■ conclusion. 

For a = 0, the system F{u,p,0) = reads: 



Pk > ei 



(24) 



giPK) = q{pk) 



VA' e M 



Since g{-) is one to one from (0, +00) to (0, +00) and thanks to the stability of the bilinear form a(-, •), this 
system has one and only one solution (which, for the pressure, reads of course pk = P*kj VAT G M.), which 
satisfies: 

\\u\\* < C3, 62 = min p*K <p < max p]^ = C4 (25) 
KeM KeM 

In addition, the jacobian of this system is block diagonal: the first block, associated to the first relation, is 
constant (this part of the system is linear) and non-singular; the second one, associated to the second relation, 
is diagonal, and each diagonal entry is equal to the derivative of gi(-), taken at the considered point. As the 
function p(-) is increasing, this jacobian matrix does not vanish for the solution of the system. Let W be defined 
by: 



W^ = {(u,p)gR^ X 



such that IIpII* u < 2max(Ci, C3) and i min(ei, 62) < p < 2max(Ci, C3)} 



The topological degree of F{-, - ,0) with respect to and W does not vanish, and by inequalities ([22]) . ([24]) . ([23]) 
and ([25]) . theorem 12.51 applies . which concludes the proof. □ 
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2.4. Some cases of application 

First of all, let us give some examples for the bilinear form a(-, •), for which the theory developped in this 
work holds. The first of them is: 



a{u,v)= / u-vdx, ||w||* = ||u||L2(n)<i , \\f\\* = ||/||L2(n)'' 

This choice for a(-, •) yields a discrete Darcy-like problem which is, up to numerical integration technicalities, the 
projection step arising in the pressure correction scheme which is considered in the present paper (see section 
[3]). Note that, in this case, the boundary condition u G Hj(rj)'^ does not make sense at the continuous level; in 
addition, the considered discretization is known to be not consistent enough to yield convergence (see remark 
13.61 hereafter) for the Darcy problem. 



The bilinear form associated to the Stokes problem provides another example of application. It may read in 
this case: 



a{u,v) ~ / Vu-Vwda;, = |u|Hi(a)'* 

Jfl.h 



or, without additional theoretical difficulties: 



a(u, w) = yu / Vu • Vu dx + ^ (V • u) (V • v) dx 



this latter form, where the real number fi > is the viscosity, corresponding to the physical shear stress tensor 
expression for a compressible flow of a constant viscosity Newtonian fluid. 

In addition, consider the steady Navier-Stokes equations, or, more generally, a time step of a (semi-)implicit 
time discretization of the unsteady Navier-Stokes equations, in which case a(-, •) and / reads: 

a{u, = T~ / pu ■ V dx + (V • pw (E) u) ■ V dx + p / Vu • Vw dx + — j (V • u) (V • v) dx 
ot Jn Jnji Jn,h 3 Jq 

where the steady case is obtained with 5t = +00, /o is the physical forcing term, p* and u* stands for known 
density and velocity fields and w is an advection field, which may be u itself (and must be u in the steady case) 
or be derived from the velocity obtained at the previous time steps. Let us suppose that the following identity 
holds: 

— / {pu — p*u*) ■ udx + / (V • (X) m) • M dx > — — / pjup — / p* 
Jn Jn,h 2dt [J 52 Jo 

which is the discrete counterpart of equation ©-(i). The algorithm considered in this paper provides an example 

where this condition is verified (see section [3]). Then the present theory applies with little modifications: in 

the proof of existence theorem 12. 6i the right hand side of the preceding equation must be multiplied by the 

homotopy parameter a (an thus this term vanishes at a = 0, which yields the problem considered in step 4 

above); the (uniform with respect to a) stability in step 1 stems from the diffusion term, and steps 2 and 3 

remain unchanged. 

Note that, in the steady state case, an additional constraint is needed for the problem to be well posed, 
namely to impose the total mass M of fiuid in the computational domain to a given value. This constraint can 
be simply enforced by solving an approximate mass balance which reads: 

r Ml 

c{h) P- +V ■ pu = 
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where \U\ stands for the measure of 17, h is the spatial discretization step and c{h) > must tend to zero with 
h, fast enough to avoid any loss of consistency. With this form of the mass balance, the theory developped here 
directly applies. 



Examining now the assumptions for the equation of state in theorem 12. 6[ we see that our results hold with 
equations of state of the form: 

g{p) = •p^l'^ or, equivalently p = p'^ , where 7 > 1 

In this case, the elastic potential is given by equation (|4]), which yields: 

P(p) = ^p'^-\ pP{p)^^p^ {=-^P) 

7— 1 7^1 7^1 

The same conclusion still holds with 7 = 1 (i.e. p = p), with P{p) = log(p) satisfying equation Q. The case 
7 > 1 is for instance encountered for isentropic perfect gas flows, whereas 7 = 1 corresponds to the isothermal 
case. It is worth noting that this range of application is larger than what is known for the continuous case, for 
which the existence of a solution is known only in the case 7 > d/2 [15, 28, 32]. 



3. A PRESSURE CORRECTION SCHEME 

In this section, we build a pressure correction numerical scheme for the solution of compressible barotropic 
Navicr-Stokcs equations ([T]), based on the low order non-conforming finite element spaces used in the previous 
section, namely the Crouzeix-Raviart or Rannachcr-Turck elements. 

The presentation is organized as follows. First, we write the scheme in the time semi-discrete setting (section 
13. ip . Then we prove a general stability estimate which applies to the discretization by a finite volume technique 
of the convection operator (section 13. 2p . The proposed scheme is built in such a way that the assumptions of 
this stability result hold (section [231); this implies first to perform a prediction of the density, as a non-standard 
first step of the algorithm and, second, to discretize the convection terms in the momentum balance equation by 
a finite volume technique specially designed to this purpose. The discretization of the projection step (section 
13. 4p also combines the finite element and finite volume methods, in such a way that the theory developped in 
section [2] applies; in particular, the proposed discretization allows to take benefit of the pressure or density 
control induced by the pressure work, i. e. to apply theorem 12.11 The remaining steps of the algorithm are 
described in section 13.51 and an overview of the scheme is given in section 13.61 The following section (section 
13. 7p is devoted to the proof of the stability of the algorithm. Finally, we shortly address some implementation 
difficulties f section [3. 8p . then provide some numerical tests (section [3. 9p performed to assess the time and space 
convergence of the scheme. 



3.1. Time semi-discrete formulation 

Let us consider a partition = to < < • • • < = of the time interval (0,T), which, for the sake of 
simplicity, we suppose uniform. Let 5t be the constant time step 5t — tk+i — tfc for = 0, 1, . . . , n. In a time 
semi-discrete setting, the scheme considered in this paper reads: 
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(26) 
(27) 



2 - Solve for : - V • (^Vp^^) = -V • {y=f^^P'') 



3 - Solve for : 

=n+l f,n+l _ „n n 

^ f- + V • (p"+i u" ® t2"+i) + Vp"+i - V • r(M"+i) = (28) 

ot 

4 - Solve for p"+\ : 

„-;n+l _ ,--,«+ 1 

— sT — + ~ ^ ° 
^(^qi^ + v. (,(,"-)."-) =0 (29) 



5 ~ Compute given by: ^p^+i ^ ^p^+i^n+i (3q) 

The first step is a prediction of the density, used for the discretization of the time derivative of the momentum. 
As remarked by Wesseling et al [2,43] , this step can be avoided when solving the Euler equations: in this case, the 
mass flowrate may be chosen as an unknown, using the explicit velocity as an advective field in the discretization 
of the convection term in the momentum balance; the velocity is then updated by dividing by the density at 
the end of the time step. For viscous flows, if the discretization of the diffusion term is chosen to be implicit, 
both the mass flowrate and the velocity appear as unknowns in the momentum balance; this seems to impeed 
the use of this trick. Let us emphasize that the way the step one is carried out (i.e. solving a discretization of 
the mass balance instead as, for instance, performing a Richardson's extrapolation) is crucial for the stability. 

Likewise, the second step is a renormalization of the pressure the interest of which is clarifled only by the 
stability analysis. A similar technique has already been introduced by Guermond and Quartapelle for variable 
density incompressible flows [18]. 

Step 3 consists in a classical semi-implicit solution of the momentum equation to obtain a predicted velocity. 

Step 4 is a nonlinear pressure correction step, which degenerates in the usual projection step as used in 
incompressible flow solvers when the density is constant (e.g. [29]). Taking the divergence of the flrst relation 
of (j29p and using the second one to eliminate the unknown velocity u"'^^ yields a non-linear elliptic problem for 
the pressure. This computation is formal in the semi-discrete formulation, but, of course, is necessarily made 
clear at the algebraic level, as described in section [3?8l Once the pressure is computed, the flrst relation yields 
the updated velocity and the third one gives the end-of-step density. 

3.2. Stability of the advection operator: a finite-volume result 

The aim of this section is to state and prove a discrete analogue to the stability identity ©-(i), which may 
be written for any sufficiently regular functions p, z and u as follows: 



/ 



^ + V.(pz.) 



zdx = — — / nz^ dx 
2 dt ' ' 



In L 

and holds provided that the following balance is satisfied by p and u: 

|+V.(p.) = 
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As stated in introduction, applying this identity to each component of the velocity yields the central argument 
of the proof of the kinetic energy theorem. 

The discrete analogue to this identity is the following. 

Theorem 3.1 (Stability of the advection operator). Let {px)KeM o-iT-d {pK)KeM be two families of positive 
real number satisfying the following set of equation: 



St 



{pk-Pk)+ E ^-,^ = 

(j=K\L 



(31) 



where F„_k is a quantity associated to the edge a and to the control volume K ; we suppose that, for any internal 
edge a = K\L, F^^k = —F„,l- Let {z*j^)k^m o.nd {zk)k<£M be two families of real numbers. For any internal 
edge a = K\L, we define Za either by = \{zk + zi), either by z^ = zk if F(j^k > and z^ ~ zl otherwise. 
The first choice is referred to as an the "centered choice" , the second one as "the upwind choice". In both cases, 
the following stability property holds: 



^ {PK ZK - p*K Zk) + E -^'^^■^ ^° 

a=K\L 



- 2 ^ 6t 



2 2 

PK Zk - Pk Zk 



(32) 



Proof. We write: 



where Ti and T2 reads: 



K<^M 



\K\ . * Mi P 

-jj- [PK ZK - Pk Zk) + 2^ ^<yJ< Za 

a=K\L 



T1+T2 



5t 



y Fa-.K Za 

r=K\L 



The first term reads: 



\K\ 



^-L= J2 ~sr (PK - Pk) + Pk ZK {zk - z*j()] 
KeM 

Developping the last term by the identity a{a — b) = ^(a^ + (a — 6)^ — 6^), we get: 

^1=2^ Zk{pk - Pk) + 1^ 2^ Pk{zk- Zk ) + 7^ 2^ PkK^k - Zk) 

KeM 



KeM 



* \2 



2 ^ 6t 

KeM 



1,3 



The last term, namely T1.3, is always non-negative and can be seen as a dissipation associated to the backward 
time discretization of equation p2p . We now turn to T2 : 



T2 = E 4 




+ E ^K 


E ^'^•K {Za - ZK 


KeM 


a=K\L 


KeM 


a=K\L 



T2A 



To 
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The first term, namely T2^i, will cancel with Ti^i by equation (PT|) . The second term reads, developping as 
previously the quantity zk {za — zk)'- 



a=K\L 



KeM 



J2 FaJ< [(^. - ^K f - zl] 

<7=K\L 



T2, 



Reordering the sum in the last term, wc have, as F^^k — ~Fa,L- 



crG£i„t (<y=K\L) 

This expression can easily be seen to vanish with the centered choice. With the upwind choice, supposing 
without loss of generality that we have chosen for the edge a = K\L the orientation such that F^^^k > 0, we 
get, as Za- = Zk- 



T, 



2.3 



ce£i„t {>y=K\L) 



We thus have, by equation (PT|) : 



T2,2 > - 



'A- 



cr=K\L 



1 V- 1^1 2 . 

2 E -w'^^p^ 



Pk) 



and thus: 



Ti+T2> 



\K\ r 



2 ^ St 



Zk (pk - Pk) + Pk {zk - z'k ) 



which concludes the proof. 



□ 



Remark 3.2. Equation ([?!]) can be seen as a discrete mass balance, with F^^k standing for the mass flux 
across the edge cr, and the right hand side of p2l) may be derived by the multiplication by zk and summation 
over the control volumes of the transport terms in a discrete balance equation for the quantity pz, reading: 



St 



(Pk zk-PkZ*j^) + 



<7=K\L 



Fa.K Za 



[possible diffusion terms] ■ 







In this context, the relation pip is known to be exactly the compatibility condition which ensures a discrete 
maximum principle for the solution z of this transport equation, provided that the upwind choice (or any 
monotone choice) is made for the expression of Za [27]. We proved here that the same compatibility condition 
ensures a stability for p^^^z. 

3.3. Spatial discretization of the density prediction and the momentum balance equation 

The main difficulty in the discretization of the momentum balance equation is to build a discrete convection 
operator which enjoys the stability property ©-(i). To this purpose, we derive for this term a finite volume 
discretization which satisfies the assumptions of theorem 13. II 

The natural space discretization for the density is the same as for the pressure, i.e. piecewise constant 
functions over each element. For the Rannacher-Turek element, this legitimates a standard mass lumping 
technique for the time derivative term, since no additional accuracy seems to have to be expected from a more 
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Figure 1. Dual finite volume mesh: the so-called "diamond cells". 

precise integration. For the Crouzeix-Raviart element, the mass matrix is already diagonal. Let the quantity 
|£>(j| be defined as follows: 

\D,\'':^ [ ip„dx>Q (33) 
Jn 

For the Crouzeix-Raviart element, \D„\ can be identified to the measure of the cone with basis a and with 
the mass center of the mesh as opposite vertex. The same property holds for the Rannacher-Turek element 
in the case of quandrangles {d = 2) or cuboids (d — i), which arc the only case considered here, even though 
extensions to non-perpendicular grids are probably possible. 

For each internal edge a — K\L, this conic volume is denoted by Dx.a', the volume = Dx.a U -Dl.o- is 
referred to as the "diamond cell" associated to cr, and DK,a is the half-diamond cell associated to a and K (see 
figure [T]). The measure of Dk,^ is denoted by \DK,a\- 

The discretization of the term p" u" thus leads, in the equations associated to the velocity on cr, to an 
expression of the form p" , where results from an average of the values taken by the density in the two 
elements adjacent to cr, weighted by the measure of the half-diamonds: 

Va e fint, \Da\ pl = \DkA Pk + \DlA pI (34) 



In order to satisfy the compatibility condition introduced in the previous section, a prediction of the density 
is first performed, by a finite volume discretization of the mass balance equation, taking the diamond cells as 
control volumes: 

Vaefint, ^(pr'-p';)+ E ^"^' = (35) 

where £{0^) is the set of the edges of D^^ and o. stands for the mass fiux across e outward Dg^. This latter 
quantity is expressed as follows: 

= |£| < • n,,, p^+i 

where |e| is the measure of e, Ug^^ is the normal to e outward Da, the velocity u" is obtained by interpolation 
at the center of e in m" (using the standard finite element expansion) and p^'^^ is the density at the edge, 
calculated by the standard upwinding technique [i.e. either p^"^^ if -n^^a > or p^t^ otherwise, with cr' such 
that £ separates D^ and Do-', which we denote by e = Da\Da>). 
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The discretization convection terms of the momentum balance equation are built from relation ([55)) , according 
to the structure which is necessary to apply theorem 13. II This yields the following discrete momentum balance 
equation: 



e(^£(D„), 
e=D„\D„, 

+ \ t(w"+1) : V(^(;) Ax- I V • Lpf dx= [ 

Jn.h Jn.h Jn 



(36) 

• V'i*^ cr e fi„t, l<i<d 



where we recall that (/s^ the vector shape function associated to Va-,i , which reads ipa '^ = (pa- with is the i*'* 
vector of the canonical basis of K.'' and ip^ the scalar shape function, and the notation ^ means 'YliKeM Ik' 

Note that, for Crouzcix-Raviart elements, a combined finite volume/finite clement method, similar to the 
technique employed here for the discretization of the momentum balance, has already been analysed for a 
transient non-linear convection-diffusion equation by Feistauer and co-workers [1, 11, 16]. 

3.4. Spatial discretization of the projection step 

The fully discrete projection step used in the described algorithm reads: 

~n+l r 

\Da\ ^ (S^r - KV) + / - V • dx = a e £i„t, l<i<d 



St 

' St 



Jn,h 



a=K\L 



where (v+^)"+i and iy'^KT""^ stands respectively for max(v""^\ 0) and max(-v^-y , 0) with v^^^^ = \cj\ u^+i • 
TiKL- The first (vector) equation may be seen as the finite element discretization of the first relation of the 
projection step (P^ . with the same lumping of the mass matrix for the Rannacher-Turek element as in the 
prediction step. As the pressure is piecewise constant, the finite element discretization of the second relation 
of (|29p . i.e. the mass balance, is equivalent to a finite volume formulation, in which we introduce the standard 
first-order upwinding. Exploiting the expression of the velocity and pressure shape functions, the first set of 
relations of this system can be alternatively written as follows: 

~n+l 

\D,\ ^ «+i - <+i) + \a\ [{pl+' ~ ipl+' -pl+')] UKL = V<T e fint, a = K\L (38) 

or, in an algebraic setting: 

^ M.,.+i (i2,"+i - {*"+!) + B* (p"+i - p"+i) = (39) 
St 

In this relation, stands for the diagonal mass matrix weighted by (wCT)ae£i„t (^Oj fo'' 1 < * < and a G £int, 
the corresponding entry on the diagonal of Mp^+i reads (Mp^+i )cr,i = (1-0^1 p^"*"^)), B* is the matrix of R^''^)^*^, 
where is the number of internal edges (i.e. N ~ card (£int)) and M is the number of control volumes in the 
mesh (i.e. M = card(A^)), associated to the gradient operator; consequently, the matrix B is associated to 
the opposite of the divergence operator. Throughout this section, we use the same notation for the discrete 
function (defined as usual in the finite element context by its expansion using the shape functions) and for 
the vector gathering the degrees of freedom; so, in relation ([39]) . u stands for the vector of M''^ components 
W(T,i, 1 < j < rf, cr G fint and p stand for the vector of R^^ of components pk, K G M. Both forms p8|) and 
are used hereafter. 



We have the following existence result. 
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Proposition 3.3. Let the equation of state g(-) be such that g(0) = 0, liniz^+oo g{z) = +oo and there exists 
an elastic potential function P(-) such that the function z l—^ zP{z) is bounded from above in (0, +oo), once 
continuously differentiable and strictly convex. Then the nonlinear system ([37]) admits at least one solution and 
any possible solution is such that p^ > 0, V-ftT Cz A4 (or, equivalently, px > 0, V-ftT £ A4). 

Proof. Let us suppose that > 0, Vcr G £int- Then the theory of section [2] apphes, with: 

This yields both the existence of a solution and the positivity of the pressure. In view of the form of the discrete 
density prediction ()35p . this latter property extends by induction to any time step of the computation (provided, 
of course, that the initial density is positive everywhere). □ 

We finish this section by some remarks concerning the projection step at hand. 
Lemma 3.4. The following identity holds for each discrete pressure q Cz L^: 

yKeM, (B M-i^, B* q)K = -^l^ilK- Ql) 



Proof. Let g G L/i be given. By relation ((38)) . we have: 

(B* q)„^i ^\(7\ [qK - qh) riKL ■ Ri 

Let Ik S Lh be the characteristic function of K. Denoting by (•, •) the standard Euclidian scalar product, by 
the previous relation and the definition of lumped velocity mass matrix, we obtain: 

(B M~ 1 B* g, Ik) = (Mri B* g, B* 1^) 

d ^ 

= E E -:K+TT7^-^\'^\'^1K-qL)nKL■e^][\(J\nKL■e^] 

aes-.^u <y=K\L i=i l^-^l 

which, remarking that X]f=i("'^i ' ^«)^ = 1; yields the result. □ 

Remark 3.5 (On spurious pressure boundary conditions). In the context of projection methods for uncom- 
pressible flow, it is known that spurious boundary conditions are to be imposed to the pressure in the projection 
step, in order to make the definition of this step of the algorithm complete. These boundary conditions are 
explicit when the process to derive the projection step is first to pose the problem at the time semi-discrete 
level and then discretize it in space; for instance, with a constant density equal to one and prescribed velocity 
boundary conditions on 917, the semi-discrete projection step would take the form: 

-A (p"+i - p"+i) = -^V • t2"+i in Vt 

ot 

V(p"+^ -p"+^) -n = on do. 

When the elliptic problem for the pressure is built at the algebraic level, the boundary conditions for the 
pressure are somehow hidden in the discrete operator BM~^B*. Lemma 13.41 shows that this matrix takes the 
form of a finite- volume Laplace discrete operator, with homogeneous Neumann boundary conditions, i.e. the 
same boundary conditions as in the time semi-discrete problem above stated. 
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Remark 3.6 (On the non-consistency of the discretization at hand for the Darcy problem). Considering the 
semi-discrctc problem p9|) . taking p^'^^ = 1, Vcr € £int, one could expect to recover a consistent discretization 
of a Poisson problem with homogeneous Neumann boundary conditions. The following example shows that this 
route is misleading. Let us take for the mesh a uniform square grid of step h. The coefficient |crp/|Dc,| can be 
easily evaluated, and we obtain: 

(BM-lB*q)K=d ^^^^K-qL) 
a=K\L 

that is the usual finite volume Laplace operator, but multiplied by the space dimension d. This result is of 
course consistent with the wcUknown non-consistency of the Rannacher-Turek element for the Darcy problem, 
and similar examples could also be given for simplicial grids, with the Crouzeix-Raviart clement. 

3.5. Renormalization steps 

The pressure renormalization (step 2 of the algorithm) is introduced for stability reasons, and its form can 
only be clarified by the stability study. It reads: 

B M-i_^i B* = B B* p" (40) 

where the density at the face at time level n necessary to the definition of M~i is obtained from the 

density field p" € Lh by equation ((M]) . In view of the expression of these operators provided by lemma [SHI this 
relation equivalently reads: 

VA^ ^M, E ^ ^ (P"/^ - PT') = E ^ - Vl) (41) 

a=K\L P'' l^'^l a=K\L VP<t ^/pi l^'^l 

As B* and B stands for respectively the discrete gradient and (opposite of the) divergence operator, this system 
can be seen as a discretization of the semi-discrete expression of step 2; note however, as detailed in remark [3.61 
that this discretization is non-consistent. 

The velocity renormalization (step 5 of the algorithm) simply reads: 



Vaefint, V/5S+i<+i or u"+' = u-+' (42) 

3.6. An overview of the algorithm 

To sum up, the algorithm considered in this section is the following one: 

(1) Prediction of the density - The density on the edges at t", {Pa)ije£int- being given by compute 
(/5"+"'^)o-G£i„t by the upwind finite volume discretization of the mass balance over the diamond cells psp . 

(2) Renormalization of the pressure - Compute a rcnormalizcd pressure {p^^)KeM by equation (|4ip . 

(3) Prediction of the velocity - Compute {u^~^^)cr^£i,^^ by equation obtained by a finite volume dis- 
cretization of the transport terms over the diamond cells and a finite element discretization of the other 
terms. 

(4) Projection step - Compute {u^a'^^)ae£-^^ and {p^^)KeM from equation ([37]) . obtained by a finite element 
discretization of the velocity correction equation and an upwind finite volume discretization of the mass 
balance (over the elements K G A4). 

(5) Renormalization of the velocity - Compute (tiJJ+^)cg£;^j from equation (|42)) . 
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The existence of a solution to step 4 is proven above; the other problems are linear, and their regularity 
follows from standard coercivity arguments, using the fact that the discrete densities {i.e. p" and are 
positive, provided that this property is satisfied by the initial condition. 

3.7. Stability analysis 

In this section, we use the following discrete norm and semi-norm: 



where p = {pa)crGeiy,t is a family of positive real numbers. The function || • \\1 - defines a norm over Wh, and 
I • |ft,_p can be seen as a weighted version of the semi- norm classical in the finite volume context [14]. The 
links between this latter semi-norm and the problem at hand are clarified in the following lemma, which is a 
straightforward consequence of lemma [ 



Lemma 3.7. The following identity holds for each discrete pressure q G Lh 

{BM^'B'q,q)^\q\l^ 



We are now in position to state the stability of the scheme under consideration. 

Theorem 3.8 (Stability of the scheme). Let the equation of state g{-) be such that g(0) = 0, lim^^^+oo p(-z) = 
-foo and there exists an elastic potential function P(-) satisfying ([3]) such that the function z t-^ z P{z) is 
bounded from above in (0, -foo), once continuously differentiate and strictly convex. Let {u"')o<n<N , iu"')o<n<N , 
(p")o<n<Af cifid {p"')o<n<N bc the solution to the considered scheme, whith a zero forcing term. Then the 
following bound holds for n < N : 



n+l 

(44) 



^iKlllpo + /^pW)dx + ^|/| 



< 

- 2 



2 

0|2 



Proof. Multiplying each equation of the step 3 of the scheme (|36p by the corresponding unknown {i.e the 
corresponding component of the velocity {t"+^ on the corresponding edge a) and summing over the edges and 
the components yields, by virtue of the stability of the discrete advection operator (theorem 13. ip : 

^ ll""+iLp"+i - ^ ll«"llip" + ^ riu-^') ■■ VS"+^ dx - j^r^'^ ■ dx < (45) 

— 1 /2 

On the other hand, reordering equation ((39)) and multiplying by M-„+i (recall that Mp„+i is diagonal), we 
obtain: 

g p"+i p" + i 1- c, P-+1 p..+ i 
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Squaring this relation gives: 
which reads: 



(Mp„+iS"+\ + (Mr„V, B*p"+\ B*p"+i) + I B*p^ 



Muhiplying by 5t/2, remarking that, Vw G VK/i, (Mp,j+iz;, u) = || i'|L ^n+i Q-iid that, thanks to lemina 13.71 
Vg e (M-^, B*g, B*(7) = (BM-^, B*g,q) = |g|^_^„+i, we get: 



^ ll'""'''^llh,p"+i + y 1/^ l/i,p"+ 



The quantity — ({t"+^, B*p"+-'^) is nothing more than the opposite of the term / p"''^^V ■ u^'^'^ dx appearing 
in (|45)l . so summing (|45)l and (|46l) makes these terms disappear, leading to: 



4 -2^11" 1"^'"" 



Finally, B*p"+^) is precisely the pressure work which can be bounded by the time derivative of the elastic 

potential, as stated in theorem |2. II 



(47) 



The proof then ends by using the renormalization steps (step 2 and 5 of the algorithm). Step 2 reads in an 
algebric setting: 

B M" B* = B MZlil Mr„i/2 B* 

Multiplying by we obtain: 

(m-Y? B*p"+i, M-^f B*|5"+i) = (Mry2B*p", M-Yf B*p"+i) 
and thus, by Cauchy-Schwartz inequality: 
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This relation yields |p"^^|^pTi+i < \p"'\'hp"- addition, step 5 of the algorithm gives p„+i = 

||it"+^||^ p"+i ■ Using these two relations in (|T7)) . we get: 

and the estimate of theorem 13.81 follows by summing over the time steps. □ 

Remark 3.9 (On the upwinding of the mass balance discretization, the inf-sup stability of the discretization 
and the appearance of spurious pressure wiggles.). In the scheme considered in this section, the upwinding in 
the discretization of mass balance controls the onset of density oscillations. As long as the pressure and the 
density are linked by an increasing function, that is as long as the flow remains compressible with a reasonable 
equation of state, it will probably be sufficient to prevent from pressure oscillations. Besides, the fourth term 
of the left hand side of (|44)) . i.e. the term involving provides a control on the discrete semi- 

norm of the pressure, at least for large time steps, and therefore also produces an additional pressure smearing. 
However, it comes up in the analysis as the composition of the discrete divergence with the discrete gradient; 
consequently, one will obtain such a smoothing effect only for inf-sup stable discretizations. Note also that, 
even for steady state problems, some authors recommend the use of stable approximation space pairs to avoid 
pressure wiggles [3,17]. 



Remark 3.10 (On a different projection step). Some authors propose a different projection step [2,43], which 
reads in time semi-discrete setting: 



5t 

It 



+ V(p"+^ -p"+^) = 



+ V • (£'(p"+^)u"+^) = 



Considering this system, one may be tempted by the following line of thought: choosing (7"+^ = £i(|)""'"^) 
as variable, taking the discrete divergence of the first equation and using the second one will cause the convec- 
tion term of the mass balance to disappear from the discrete elliptic problem for the pressure, whatever the 
discretization of this term may be. Consequently, the equation for the pressure will be free of the non-linearities 
induced by the upwinding and the dependency of the convected density on the pressure, while one still may 
hope to obtain a positive upwind (with respect to the density) scheme. In fact, this last point is incorrect. To 
be valid, it would necessitate that, from any solution ((7"+^,^"+^), one be able to compute a velocity field 
by dividing g"^^ by the density of the control volume located upstream with respect to u""*"^. Unfortunately, 
it is not always possible to obtain this upstream value; for instance, if for two neighbouring control volumes K 
and L, pK < 0, > and g""*"^ • nK\L > 0, neither the choice of K nor L for the upstream control volume 
is valid. Consequently, with this discretization, we are no longer able to guarantee neither the positivity of p 
nor the absence of oscillations. However, as explained below, if the density remains positive, we will have a 
smearing of pressure or density wiggles due to the fact that the discretization is inf-sup stable. 

3.8. Implementation 

The implementation of the first three steps P51) - ([^51) of the numerical scheme is standard, and we therefore 
only describe here in details the fourth step, that is the projection step. The precise algebraic formulation of 
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the system ([^ reads: 

1 

6i 



Mp„+i (t2"+i - + B* - = 



(48) 



where M^n+i and Q^n+i are two diagonal matrices; for the first one, we recall that the entry corresponding to 
an edge a £ £int, cr = K\L is computed by multiplying the measure of the diamond associated to a by the 
predicted density (at the edge center) p^"*"^; in the second one, the same entry is obtained by just taking the 
density at in the element located upstream of a with respect to i.e. either g{p^^) or g{p^^^). 

Note that these definitions can be extended in a straightforward way for the boundary edges, if the velocity is 
not prescribed to zero on the boundary of the computational domain. The matrix R is diagonal and, for any 
K E M, its entry Kk is the measure of the element K. For the sake of simplicity, we suppose for the moment 
that the equation of state is linear: 



The elliptic problem for the pressure is obtained by multiplying the first relation of (jlS]) by B Q^n+i (Mp^ 
and using the second one. This equation reads: 



where L = BQ n+i (M^n+i)^^ B* can be viewed, for the discretization at hand, as a finite volume discrete 

Pup ^ r 

approximation of the Laplace operator with Neumann boundary conditions (when the velocity is prescribed at 
the boundary), multiplied by the space dimension d and the densities ratio (see remarks 13. 51 and 13. 6p . We recall 
that, by a calculation similar to the proof of lemma [3^ this matrix can be evaluated directly in the "finite 
volume way" , by the following relation, valid for each element K: 

where Pup.o- stands for the upwind density associated to the edge cr. Provided that p"'''^ is known, the first 
relation of gives us the updated value of the velocity: 

^"+1 ^ ^n+l _ (Mp„ + l)-l B* (p"+l - p"+^) 

As, to preserve the positivity of the density, we want to use in the mass balance the value of the density 
upwinded with respect to i2"+^, equations p9|) and (|3.8p are not decoupled, by contrast with what happens in 
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usual projection methods. We thus implement the following iterative algorithm: 
Initiahzation: p^^^ ~ and Wg 



Step 4.1 - Solve for p^+l 



dp St^ 



/2 
n+l 



Pk+1/2 - 



where the density in L and Q^^+i is evaluated at p? and the upwinding 

pup 

in Q„n+i is performed with respect to 

Pup 



n+l 



Step 4.2 - Compute p^+l as p'^+l = ap^+l^^ + (1 - a) pl+^ 

Step 4.3 - Compute u)!^^ as : 

-5t (M,„+0-'B*(p^+i-p"+i) 



Convergence criteria: max [ Wp^'^l — P^^^ || , 



.'1+1 



< £ 



The second step of the previous algorithm is a relaxation step which can be performed to ensure convergence; 
however, in the tests presented hereafter, we use a = 1 and obtain convergence in few iterations (typically less 
than 5). When the equation of state is nonlinear, step 4.1 is replaced by one iteration of Newton's algorithm. 

3.9. Numerical experiments 

In this section, we describe numerical experiments performed to assess the behaviour of the pressure cor- 
rection scheme presented in this paper, in particular the convergence rate with respect to the space and time 
discretizations. 

With Q = (0, 1) X (— ^, \), we choose for the momentum and density the following expressions: 



pu = cos(7ri) 



sin(7rx^"'"^) 
cos(7ra;'^^^) 



p — 1 + - sin(7ri) cos(7ra;'-"'"-') — sin(7ra;'-^'') 

These functions satisfy the mass balance equation; for the momentum balance, we add the corresponding right 
hand side. In this latter equation, the divergence of the stress tensor is given by: 

V ■ t{u) = pAu+ -W ■ u, p. = 10^'^ 
3 

The pressure is linked to the density by the following equation of state: 

p-1 



7 = 1.4, Ma = 0.5 



7 Ma^ 

where the parameter Ma corresponds to the characteristic Mach number. 

We use in these tests a special numerical integration of the forcing term of the momentum balance, which 
is designed to ensure that the discretization of a gradient is indeed a discrete gradient {i.e. if the forcing term 
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Figure 2. velocity error as a function of the time step. 




Figure 3. pressure error as a function of the time step. 



/ can be recast under the form / = Vg, the discrete right hand side of the momentum balance belongs to the 
range of B'). 

Velocity and pressure errors obtained at t = 0.5, in respectively and discrete norms and as a function 
of the time step, are drawn on respectively figure[2]and figure[31 for 20 x 20, 40 x 40 and 80 x 80 uniform meshes. 
For large time steps, these curves show a decrease corresponding to approximately a first order convergence in 
time for the velocity and the pressure, until a plateau is reached, due to the fact that errors are bounded from 
below by the residual spatial discretization error. For both velocity and pressure, the value of the errors on this 
plateau show a spatial convergence order (in norm) close to 2. 

4. Conclusion 

We presented in this paper a numerical scheme for the barotropic Navier-Stokes compressible equations, 
based on a pressure-correction time stepping algorithm. For the spatial discretization, it combines low-order 
non-conforming mixed finite elements with finite volumes; in the incompressible limit, one recovers a classical 
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projection scheme based on an inf-sup stable pair of approximation spaces for the vefocity and the pressure. 
This scheme is proven to enjoy an unconditional stability property: irrespectively of the time step, the discrete 
solution obeys the a priori estimates associated to the continuous problem, i. e. strict positivity of the density, 
bounds in i°°-in-timc norm of the quantity pu^ dx and J^^ pP{p) dx and in i^-in-timc norm of the viscous 
dissipation t(u) : Vwdx. To our knowledge, this result is the first one of this type for barotropic compressible 
flows. 

However, the scheme presented here is by no means "the ultimate scheme" for the solution to the compressible 
Navier-Stokes equations. It should rather be seen as an example of application (and probably one of the less 
sophisticated ones) of the mathematical arguments developped to obtain stability, namely theorems 12 . II (discrete 
elastic potential identity) and l3.1l fstabilitv of the advection operator), and our hope is that these two ingredients 
could be used as such or adapted in the future to study other algorithms. For instance, a computation close 
to the proof of theorem 13.81 (and even simpler) would yield the stability of the fully implicit scheme; adding to 
this latter algorithm a prediction step for the density (as performed here) would also allow to linearize (once 
again as performed here) the convection operator without loss of stability. A stable pressure-correction scheme 
avoiding this prediction step can also be obtained, and is currently under tests at IRSN for the computation 
of compressible bubbly flows. Besides these variants, less diffusive schemes should certainly be searched for. 
Finally, the proposed scheme is currently the object of more in-depth numerical studies including, in particular, 
problems admitting less smooth solutions than the test presented here. 

References 

[1] Ph. Angot, V. Dolejsf, M. Feistauer, and J. Felcman. Analysis of a combined barycentric finite volume-nonconforming finite 
element method for nonlinear convection-diffusion problems. Applications of Mathematics, 4:263-310, 1998. 

[2] H. Bijl and P. Wesseling. A unified method for computing incompressible and compressible flows in boundary-fitted coordinates. 
Journal of Computational Physics, 141:153-173, 1998. 

[3] M.O. Bristeau, R. Glowinski, L. Dutto, J. Periaux, and G. Roge. Compressible viscous flow calculations using compatible finite 
element approximations. International Journal for Numerical Methods in Fluids, 11:719-749, 1990. 

[4] V. CasuUi and D. Greenspan. Pressure method for the numerical solution of transient, compressible fluid flows. International 
Journal for Numerical Methods in Fluids, 4:1001-1012, 1984. 

[5] A.J. Chorin. Numerical solution of the Navier-Stokes equations. Mathematics of Computation, 22:745-762, 1968. 

[6] P. G. Ciarlet. Handbook of numerical analysis volume II : Finite elements methods - Basic error estimates for elliptic problems. 
In P. Ciarlet and J.L. Lions, editors, Handbook of Numerical Analysis, Volume II, pages 17—351. North Holland, 1991. 

[7] P. Colella and K. Pao. A projection method for low speed flows. Journal of Computational Physics, 149:245-269, 1999. 

[8] M. Crouzeix and P.-A. Raviart. Conforming and nonconforming finite element methods for solving the stationary Stokes 
equations I. Revue Frangaise d'Automatique, Informatique et Recherche Operationnelle (R.A.I.R.O.), R-3:33-75, 1973. 

[9] K. Deimling. Nonlinear Functional Analysis. Springer, New- York, 1980. 
[10] I. Demirdzic, Z. Lilek, and M. Peric. A collocated finite volume method for predicting fiows at all speed. International Journal 

for Numerical Methods in Fluids, 16:1029-1050, 1993. 
[11] V. Dolejsf, M. Feistauer, J. Felcman, and A. Klikova. Error estimates for barycentric finite volumes combined with noncon- 
forming finite elements applied to nonlinear convection-diffusion problems. Applications of Mathematics, 47:301-340, 2002. 
[12] A. Ern. Aide Memoire Elements finis. Dunod, Paris, 2005. 

[13] R. Eymard, T. Gallouet, M. Ghilani, and R. Herbin. Error estimates for the approximate solutions of a nonlinear hyperbolic 
equation given by finite volume schemes. IMA Journal of Numerical Analysis, 18(4):563-594, 1998. 

[14] R. Eymard, T Gallouet, and R. Herbin. Finite volume methods. In P. Ciarlet and J.L. Lions, editors, Handbook of Numerical 
Analysis, Volume VII, pages 713-1020. North Holland, 2000. 

[15] E. Feireisl. Dynamics of viscous compressible flows, volume 26 of Oxford Lecture Series in Mathematics and its Applications. 
Oxford University Press, 2004. 

[16] H. Feistauer, J. Felcman, and I. Straskraba. Mathematical and computational methods for compressible fiows. Oxford Science 

Publications. Clarendon Press, 2003. 
[17] M. Fortin, H. Manouzi, and A. Soulaimani. On finite element approximation and stabilization methods for compressible viscous 

fiows. International Journal for Numerical Methods in Fluids, 17:477-499, 1993. 
[18] J.-L. Guermond and L. Quartapelle. A projection FEM for variable density incompressible fiows. Journal of Computational 

Physics, 165:167-188, 2000. 

[19] J.L. Guermond, P. Minev, and J. Shen. An overview of projection methods for incompressible flows. Computer Methods in 
Applied Mechanics and Engineering, 195:6011—6045, 2006. 



28 TITLE WILL BE SET BY THE PUBLISHER 

[20] F.H. Harlow and A. A. Amsdcn. Numerical calculation of almost incompressible flow. Journal of Computational Physics, 
3:80-93, 1968. 

[21] F.H. Harlow and A. A. Amsden. A numerical fluid dynamics calculation method for all flow speeds. Journal of Computational 
Physics, 8:197-213, 1971. 

[22] R.I. Issa. Solution of the implicitly discretised fluid flow equations by operator splitting. Journal of Computational Physics, 
62:40-65, 1985. 

[23] R.I. Issa, A.D. Gosman, and A. P. Watkins. The computation of compressible and incompressible recirculating flows by a 

non-iterative implicit scheme. Journal of Computational Physics, 62:66-82, 1986. 
[24] R.I. Issa and M.H. Javarcshkian. Pressure-based compressible calculation method utilizing total variation diminishing schemes. 

AIAA Journal, 36:1652-1657, 1998. 
[25] K.C. Karki and S.V. Patankar. Pressure based calculation procedure for viscous flows at all speeds in arbitrary configurations. 

AIAA Journal, 27:1167-1174, 1989. 
[26] M.H. Kobayashi and J.C.F. Pereira. Characteristic-based pressure correction at all speeds. AIAA Journal, 34:272-280, 1996. 
[27] B. Larrouturou. How to preserve the mass fractions positivity when computing compressible multi-component flows. Journal 

of Computational Physics, 95:59-84, 1991. 
[28] P.-L. Lions. Mathematical topics in fluid mechanics — volume 2 - compressible models, volume 10 of Oxford Lecture Series in 

Mathematics and its Applications. Oxford University Press, 1998. 
[29] M. Marion and R. Temam. Navier-Stokes equations: Theory and approximation. In P. Ciarlet and J.L. Lions, editors. Handbook 

of Numerical Analysis, Volume VI. North Holland, 1998. 
[30] F. Moukalled and M. Darwish. A high-resolution pressure-based algorithm for fluid flow at all speeds. Journal of Computational 

Physics, 168:101-133, 2001. 

[31] P. Nithiarasu, R. Codina, and O.C. Zienkiewicz. The Characteristic-Based Split (CBS) scheme - a unifled approach to fluid 

dynamics. International Journal for Numerical Methods in Engineering, 66:1514—1546, 2006. 
[32] A. Novotny and I. Straskraba. Introduction to the mathematical theory of compressible flow, volume 27 of Oxford Lecture 

Series in Mathematics and its Applications. Oxford University Press, 2004. 
[33] G. Patnaik, R.H. Guirguis, J. P. Boris, and E.S. Oran. A barely implicit correction for flux-corrected transport. Journal of 

Computational Physics, 71:1-20, 1987. 
[34] E.S. Politis and K.C. Giannakoglou. A pressure-based algorithm for high-speed turbomachinery flows. International Journal 

for Numerical Methods in Fluids, 25:63-80, 1997. 
[35] R. Rannacher and S. Turek. Simple nonconforming quadrilateral Stokes element. Numerical Methods for Partial Differential 

Equations, 8:97-111, 1992. 

[36] R. Temam. Sur I'approximation de la solution des equations de Navier-Stokes par la methode des pas fractionnaircs II. Arch. 

Rat. Mech. Anal, 33:377-385, 1969. 
[37] D.R. van der Heul, C. Vuik, and P. Wesscling. Stability analysis of segregated solution methods for compressible flow. Applied 

Numerical Mathematics, 38:257-274, 2001. 
[38] D.R. van der Heul, C. Vuik, and P. Wesseling. A conservative pressure-correction method for flow at all speeds. Computer & 

Fluids, 32:1113-1132, 2003. 

[39] J. P. Van Dormaal, G.D. Raithby, and B.H. McDonald. The segregated approach to predicting viscous compressible fluid flows. 

Transactions of the ASME, 109:268-277, 1987. 
[40] D. Vidovic, A. Segal, and P. Wesseling. A superlinearly convergent Mach-uniform finite volume method for the Euler equations 

on staggered unstructured grids. Journal of Computational Physics, 217:277-294, 2006. 
[41] C. Wall, CD. Pierce, and P. Moin. A semi-implicit method for resolution of acoustic waves in low Mach number flows. Journal 

of Computational Physics, 181:545-563, 2002. 
[42] I. Wenneker, A. Segal, and P. Wesseling. A Mach-uniform unstructured staggered grid method. International Journal for 

Numerical Methods in Fluids, 40:1209-1235, 2002. 
[43] P. Wesseling. Principles of computational fluid dynamics, volume 29 of Springer Series in Computational Mathematics. 

Springer, 2001. 

[44] O.C. Zienkiewicz and R. Codina. A general algorithm for compressible and incompressible flow - Part I. The split characteristic- 
based scheme. International Journal for Numerical Methods in Fluids, 20:869—885, 1995. 



